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Abstract 

A complete analysis is presented for the far-field creeping flow produced by a 
multipolar force distribution in a fluid confined between two parallel planar walls. 
We show that at distances larger than several wall separations the flow field as- 
sumes the Hele-Shaw form, i.e., it is parallel to the walls and varies quadratically 
in the transverse direction. The associated pressure field is a two-dimensional har- 
monic function that is characterized by the same multipolar number m as the orig- 
inal force multipole. Using these results we derive asymptotic expressions for the 
Green's matrix that represents Stokes flow in the wall-bounded fluid in terms of 
a multipolar spherical basis. This Green's matrix plays a central role in our re- 
cently proposed algorithm [Physica A xx, xxx (2005)] for evaluating many-body 
hydrodynamic interactions in a suspension of spherical particles in the parallel-wall 
geometry. Implementation of our asymptotic expressions in this algorithm increases 
its efficiency substantially because the numerically expensive evaluation of the exact 
matrix elements is needed only for the neighboring particles. Our asymptotic anal- 
ysis will also be useful in developing hydrodynamic algorithms for wall-bounded 
periodic systems and implementing acceleration methods by using corresponding 
results for the two-dimensional scalar potential. 

Key words: hydrodynamic interactions, confined systems, Stokes flow, 
suspensions, Hele-Shaw flow 



1 Introduction 



Numerical and theoretical investigations of particle motion in suspensions 
bounded by planar walls require efficient methods for evaluating hydrody- 
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namic interactions in these systems. Examples of phenomena where the hy- 
drodynamic wall effects are important include collective particle motion in 
quasi-bidimensional colloidal suspensions [3, Q, El, 0, El , and conformation dy- 
namics of a DNA molecule in a parallel-plate microchannel 0. 

Several methods for evaluating hydrodynamic interactions in wall-bounded 
systems have been proposed. In some studies, the flow reflected from the 
walls was calculated numerically using either boundary-integral 0] or finite- 
difference techniques. In a different approach (if, the exact point-force 
solution for the flow between the walls [9( was used. Wall effects were also 
included using a multiple- reflection technique [nj , and several approximation 
methods were proposed 11, 0, 0|. While all of these methods have their 



merits, they also have some essential disadvantages, such as a high numerical 
cost or an insufficient (in many cases unknown) accuracy. 

Recently we have derived 14, a novel algorithm for evaluating hydrody- 
namic friction matrix in a wall-bounded suspension of spheres under creeping- 
flow conditions. 3 Our Cartesian-representation method relies on transforma- 
tions between a spherical basis set of solutions of Stokes equations (this set 
is consistent with the particle geometry) and a Cartesian basis set (which is 
consistent with the wall geometry). The algorithm provides highly accurate 
results for multiparticle friction and mobility matrices. 

Using our approach, we have obtained several interesting numerical results. 
In particular, we have shown that the friction matrix undergoes a crossover 
from the quasi-three-dimensional to quasi-two-dimensional form when the in- 
terparticle distance becomes larger than the wall separation H. We have also 
observed an unusually large resistance coefficient for a long rigid chain of 
spheres in transverse motion (with respect to the orientation of the chain) 
in a narrow, wall-bounded space. Since both these effects involve flow on the 
length scale I ^> H, they are not captured by the usual single- wall super- 
position approximation which does not properly take the far-field flow into 
account (as demonstrated in lli 



Large-scale studies of particle dynamics in the two-wall geometry require ef- 
ficient simulation algorithms. In our approach [lil . Il5l ] the most expensive 
part is evaluation of the Green's matrix G in the multipolar representation. 
This matrix is a key quantity in our algorithm — its elements correspond to 
the coefficients in the expansion of the hydrodynamic Green's tensor for the 
wall-bounded system into multipolar basis fields. The inverse of the Green's 
matrix combined with the one-particle reflection matrices yields the multi- 
particle hydrodynamic friction matrix. 



In our algorithm [1J, |Lj| the matrix G is expressed in terms of lateral Fourier 
3 An algorithm based on similar ideas was also developed by Jones [1(J. 
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integrals with respect to the two-dimensional wave vector in a plane parallel 
to the walls. Evaluation of these integrals is especially difficult for widely 
separated particles due to the oscillatory character of the integrands. In the 
present paper we derive much simpler asymptotic formulas for the matrix G 
in the far-field regime. When the particle separation is sufficiently large, these 
formulas can be used instead of the Fourier integrals, resulting in a significant 
reduction of numerical cost (and in other important simplifications). 

Our analysis of the asymptotic form of the matrix G relies on the observa- 
tion that in the far-field regime the velocity field in the space between the 
walls assumes a simple Hele-Shaw (i.e. the lubrication) form. Accordingly, the 
flow field has only the lateral components and it varies quadratically across 
the space between the walls. Such a flow field is entirely determined by the 
corresponding pressure field, which is a two-dimensional harmonic function 
that depends only on the lateral coordinates. It follows that at large distances 
r ^> H, the full three-dimensional hydrodynamic problem is reduced to a 
much simpler two-dimensional scalar problem for the pressure. 

This paper is organized as follows. Our method 0, 0] for evaluating many- 
particle hydrodynamic interactions in the parallel-wall geometry is summa- 
rized in Sees. 2 and 3. Section 2 recalls the induced-force formulation of the 
problem, and Sec. 3 summarizes the force-multipole expansion method. The 
main theoretical results of the present analysis are given in Sec. 4, where 
the Hele-Shaw approximation for the far-field flow is discussed, and explicit 
expressions for Green's matrix G are derived. In Sec. 5 we present some re- 
sults of numerical calculations. We show the dependence of Green's matrix 
elements on the interparticle distance, and we illustrate the role of their far- 
field behavior in the description of hydrodynamic interactions in rigid arrays 
of spheres. Concluding remarks are given in Sec. 6, and some technical details 
are presented in the appendices. 



2 Multiparticle hydrodynamic interactions 

2.1 Hydrodynamic resistance 

We consider the motion of N spherical particles of the radius a, which are 
suspended in a fluid of viscosity rj, under creeping-flow conditions. The system 
is bounded by two planar parallel walls at the positions z = and z = H, 
where r = (x, y, z) are the Cartesian coordinates. The centers of particles i = 
1 . . . N are at positions Rj = (Xj, Y{, Zj), and the translational and rotational 
particle velocities are U$ and Q, j. The external forces and torques acting on the 
particles are denoted by jFj and TV It is assumed that the flow field satisfies 
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the no-slip boundary conditions on the particle surfaces and the walls. 



For a system of spheres undergoing translational and rotational rigid-body 
motion with no external flow, the particle dynamics is characterized by the 
resistance matrix 



/■tt /-tr 



Crt /rr 
ij Sjj 



hJ 



l,...,iV, 



defined by the linear relation 



N 

£ 



/"-tt ftr 

Sij Sij 

Crt /-tt 

ij Sij 



(2) 



The dot in the above equation denotes the matrix multiplication and contrac- 
tion of the Cartesian tensorial components of the resistance matrix. Our goal 
is to calculate the resistance matrix £, or its inverse, the mobility matrix fi. 
Our method 14 , 1^ for evaluating these quantities is outlined below. 



2.2 Induced-force formulation 



The effect of the suspended particles on the surrounding fluid can be described 
in terms of the induced-force distributions on the particle surfaces. These 
distributions can be written in a form 



F i (r)=5 a s ( r _ Rj )f j(r ), 



where 



5 s a (r)=a- 2 S(r-a). 
By the definition of the induced force [13, EH E^j , the flow field 



v r 



N 

£/T(r,r').F 4 (r0dr' 

i=l 



(3) 
(4) 

(5) 



is identical to the velocity field in the presence of the moving particles. Here 

T(r,r') = T (r-r / ) + T'(r,r / ) ( 6 ) 

is the Green's function for the Stokes flow in the presence of the boundaries; 
the Green's function T(r, r') is decomposed into the Oseen tensor 



Tnfr) 



(7) 
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and the part T'(r, r') that describes the flow reflected from the walls. In Eq. 
(5) it is assumed that the particles move with given velocities, but no external 
flow is imposed. 

The resistance relation (2) is linked to the induced-force distributions (3) 
through the expressions 

Ti = j F<(r) dr, Ti = jn X Fi{r) dr (8) 

for the total force and torque, respectively. To determine the resistance matrix 
(1) we thus need to evaluate the induced forces (3) for given translational and 
angular velocities of the particles. 



2. 3 Boundary-integral equations for the induced forces 



For a system of particles moving with the translational and angular veloci- 
ties Uj and f2j, the induced-force distribution (3) can be obtained from the 
boundary-integral equation 



[Z^Kr) + £ / [<%T'(r - r') + (1 - <y T(r - r')] • F,(r') dr' = vf (r), 

r E Si, (9) 

where 

< b (r) = U< + n i Xr t (10) 

is the rigid-body velocity field associated with the particle motion, and Si is 
the surface of particle i. In the boundary-integral equation (9), Zj denotes the 
one-particle scattering operator that describes the response of an individual 
particle to an external flow in an unbounded space. This operator is defined 
by the linear relation 

F^-Z^-vf), (11) 



where v™ is the velocity incident to particle i. For specific particle models (e.g., 
rigid particles or drops), explicit expressions for the operator Zj are known 
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3 Force-multipole expansion 



3. 1 Spherical basis fields 



As in a standard force-multipole approach |23j,|24j the boundary-integral equa- 
tion (9) is transformed into a linear matrix equation by projecting it onto a 
spherical basis of Stokes flow. To this end we use the reciprocal basis sets 
defined by Cichocki et al. |21]. We introduce, however, a slightly different 
normalization to exploit the full symmetry of the problem. 

The singular and nonsingular spherical basis solutions of Stokes equations 
v zma0) and V /Lx( r ) ( with / = 1, 2, . . .; m = . . . , /; and a = 0, 1, 2) have 
the following separable form in the spherical coordinates r = (r, 9, 4>): 



liri: 



CT (r)=Vr m(T (MK (/+CT) , (12a) 



Ima \ 



V 



J+a-l 



Ima ' 



The coefficients V 



Ima \ 



and V z + 



(12b) 

are combinations of vector spherical 



harmonics with angular order I and azimuthal order m. This property and the 



r-dependence in equations (12) define the Stokes- flow fields v w ( 
normalization. 



up to 



We use here a convenient normalization introduced in which emphasizes 
various symmetries of the problem. Explicit expressions for the functions V im<J 
in this normalization are given in Appendix A. We note that both in our 



present and in the original normalization 
identity HT 



21 



the basis fields v /mcr satisfy the 



^To(r-r') 



H V /m CT ( r ) V /m* CT ( r ')> 



Imo 



E v w( r ) v w( r " 



k Ima 



ty ~~~~~ :r rjrt 



13a, b) 



where T (r — r') is the Oseen tensor (7). Relation (13) assures that the Lorentz 
reciprocal symmetry of Stokes flow is reflected in the symmetry of the resulting 
matrix representation of the problem 



Following 21 1 we also introduce the reciprocal basis fields w^ cr (r), defined by 
the orthogonality relations of the form 



($t W hna I V !W) - Sll'$mm'5aa>- 



(14) 
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Here 

(A | B) = J A*(r) -B(r)dr (15) 

is the inner product, the asterisk denotes the complex conjugate, and <5„ is 
defined in Eq. (4). The reciprocal basis fields and the bra-ket notation (15) 
allows us to conveniently represent expansions of Stokes flow fields into the 
complete sets of nonsingular and singular basis fields (12). In particular, any 
Stokes flow u(r) that is nonsingular in the neighborhood of a point r = Rj 
has an expansion 

u(r) = £ v+ mCT (r - R,)<5 a s «w+ mCT « I u), (16) 

Ima 

where 

wLtW = Wma( r ~ Ri), (17a) 

^(O^^r-H,). (17b) 

3.2 Matrix representation 

The matrix representation of the boundary-integral equation (9) is obtained 
using the multipolar expansion 

F 4 (r) = £ mma)5l{v - R,)w+ m(T (r - R 4 ) (18) 

Ima 

of the induced-force distributions (3). The multipolar moments in the above 
expression are given by the projection 

Mima) = (v+ mCT W | Fj), (19) 

according to the orthogonality condition (14). The definition (19) of the mul- 
tipolar expansion is justified by the identity 

VwM = V J T (r - r')5 a S (r')w+ tCT (r') dr', (20) 

which follows from the representation (13) of the Oseen tensor. Equations (18) 
and (20) indicate that the multipolar moments fiilma) are identical (apart 
from the trivial factor rj) to the expansion coefficients of the flow field scattered 
by an isolated particle in unbounded space into the singular basis fields v^ CT . 

To obtain a linear matrix equation for the set of force multipolar moments 
fi(lma), representation (18) is inserted into the boundary-integral equation 
(9), and the resulting expression is expanded into the nonsingular basis fields 
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(12b), which yields 

N 

£ Molina | I'm'a') ft (I'm' a') = Ci {lma). (21) 

j=l I'm'a' 

For a particle moving in a quiescent fluid, the coefficients 

Ci (lma) = <5 a s «w+ a (i) | vf ) (22) 

on the right side are nonzero only for I = 1 and a = 0, 1. The matrix M in 
Eq. (21) consists of three contributions corresponding to the three terms on 
the left side of Eq. (9), 

Mij{lma | I'm'a') = (%5;/'<5 mm ' Z^~ (I; a \ a') + SyG'^tyma \ I'm'a') 

+ (1 - S^G^lma | I'm'a'). (23) 

Using the bra-ket notation these contributions can be expressed in the form 

Z-\l; a | a') = (<5 s (z)w+ CT (z) | Z- 1 | ^w^)), (24) 

G'Sma | I'm'a') = (5 s a (z)w+ m ^) | T' | 5 a s (z)w+ mV ,(z)), (25) 

and 

G^lma | I'm'a') = <5 a s «w+ mCT « I T | 6 s a (j)w+ mla ,(j)). (26) 
The first term (I; a \ a') is associated with the one particle operator Z" 1 in 
equation (9), and it relates the force multipoles fiil'm'a') induced on particle 
i to the coefficients in the expansion of the flow field incoming to this particle 
into the nonsingular spherical basis fields (12b). By the spherical symmetry, 
this term is diagonal in the indices / and m and is independent of m. The 
Green's matrices G'^ilma \ I'm'a') and Gijilma \ I'm'a') are associated with 
the integral operators that involve the kernels T'(r, r') and T(r, r'). Using the 
orthogonality relations (14) one can show that the elements of these matrices 
correspond to the expansion of the flow produced by a force multipole centered 
at Rj into the nonsingular basis (12b) centered at Rj. 

Explicit expressions for the single-particle reflection matrix Z~ l are well known 
|2lL l26| . Quadrature formulas for the Green's matrix G^ have been derived 
in our recent publication where the matrix elements Gijilma \ I'm'a) 
are represented as a combination of the free-space Green's matrix [2^, E4j 
and the wall contribution G'^ (Ima \ I'm'a) that is given in a form of a Hankel 
transform of a product of several simple matrices. The Hankel transform arises 
from angular integration of lateral Fourier modes of Stokes flow. 

The many-particle resistance matrix (1) can be obtained by solving Eq. (21) 
and projecting the induced force multipoles onto the total force and torque 
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(8). Explicit expressions for the resistance matrix in terms of the generalized 
friction matrix M~ l are given in ^^|. In numerical applications, the system of 
linear equations (21) is truncated at a given multipolar order I, and the result- 
ing approximate friction matrix is supplemented with a lubrication correction 
(as described in (l5j|). 



4 Far- field approximation 

Calculation of the exact matrix elements Gij(lma \ I'm' a') by our Cartesian- 
representation method requires numerical evaluation of Hankel transforms 
that involve the Bessel functions J m - m '{kQij)- Here k is the magnitude of the 
lateral wave vector, and Qij = \q% — Qj\, where Qi = (Xi,Yi) denotes the 
lateral position of particle i. For large interparticle distances g^, the factor 
Jm-m'(kQij) undergoes rapid oscillations as a function of k. Thus, evaluation 
of the Fourier integrals in the Hankel transforms is numerically expensive for 
such configurations. 

In the following sections we derive explicit asymptotic expressions for the 
matrix elements Gij(lma \ I'm' a') at large interparticle distances 3> H. 
As we will see, these expressions have a very simple form, and do not require 
evaluation of the Fourier integrals. 

4- 1 Hele-Shaw form of the far-field flow 

Our asymptotic analysis relies on the observation that in the far-field regime 
the flow between two parallel walls assumes the Hele-Shaw form. Accordingly, 
the asymptotic pressure field p = p as varies only in the lateral direction, and 
the associated flow field has the lubrication form 

u as (r) = -\r,- l z{H - z)V\\p as (p), (27) 

where 

r = p + ze z (28) 

and V|| is the two-dimensional gradient operator with respect to the lateral 
position p = (x, y). By the incompressibility of the flow field (27), the pressure 
field p as satisfies the two-dimensional Laplace's equation 

Vjjp as (p) = 0. (29) 

The asymptotic expressions (27) and (29) can be obtained j2?| by expanding 
the boundary-value problem for Stokes flow in the space between the walls in 
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the small parameter H/p <C 1, where p is the distance from the force distribu- 
tion that generates the fluid motion. Since the velocity field (27) itself satisfies 
the Stokes equations and boundary conditions exactly, one can show that the 
higher-order terms in the asymptotic expansion vanish. This property indi- 
cates that the correction terms are subdominant [28], which in turn suggests 
that the asymptotic behavior (27) and (29) is approached exponentially. This 
conclusion is consistent with the direct analysis of the asymptotic form of the 
Green's function in the space between the walls by Liron and Mochon [9j (see 
the discussion in Sec. 4.4 below). 



4-2 Asymptotic basis sets 



To find the far-field form of the velocity field produced by induced-force multi- 
poles and to obtain the corresponding asymptotic expressions for the elements 
of the Green's matrix Gij(lma \ I'm' a'), it is convenient to define appropri- 
ate basis sets of Hele-Shaw flow and pressure fields. The sets of singular and 
nonsingular pressures are defined by the relation 

qZHp)=V<{p), (30) 
where m = 0, ±1, ±2, . . ., and 

$ -(p) = -lnp, 3- ( p ) = -i-^-M e?»* m^O, (31a) 



$+(p) = p H e ^ (31b) 

are the two-dimensional harmonic basis functions. The associated Hele-Shaw 
basis velocity fields are 



v- ± (r) = -Hi/-^)V||$±(p), 



(32) 



according to Eq. (27). 



Below we list several useful relations for the harmonic functions (31). First, 
we have the diagonal representation for the Green's function 



J2 $-(p)$+V)> p>p'> 



m=— oo 



(33a, b) 



$+(p)<&-*(p'), p <p', 

n=— oo 

which is analogous to the representation (13) of the Oseen tensor. Next, we 
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also have the displacement theorem 

oo 

KAP ~ Qj) = E - Qi)StyT{Qa\ m I rn% (34) 

m=~ oo 

where Qij = Qi — Qj, and the displacement matrix is given by 



S+i(e-m I m') = g(- mm 0(-l)^ M±^ $-,_ m (g). (35) 



We note that due to the presence of the Heaviside step function 

9{x) -- 



0, x < 0, 

1, £>0, 



(36) 



in Eq. (35), the scalar fields with the same sign of the indices m, m' ^ do 
not couple in the displacement relation (34). We also note that the matrix 
(35) satisfies the symmetry relation 

S+i(g;m\ rm!) = S£*(-Q;m' \ m). (37) 



As a direct consequence of the displacement theorem (34) for the scalar pres- 
sure fields, we have the corresponding displacement relation for the Hele-Shaw 
basis flows (32) 

oo 

v-r(r - Qj ) = £ ' <+(r - Qi )S^{Qij\ m | mf). (38) 



m=— oo 



The term with m = in the above relation vanishes because Vq = ac- 
cording to Eqs. (31b) and (32). The prime at the summation sign has been 
introduced to emphasize that this term is omitted. 

In the following section we will derive a diagonal representation (analogous to 
(13) and (33)) for the hydrodynamic Green's tensor describing the asymptotic 
far-field response of the fluid confined between walls to a point force. 



4-3 Asymptotic Green's tensor 



An explicit expression for the far-field flow produced by a point force in the 
space between the walls has been derived by Liron and Mochon (see also 
[29|). According to their results, the far-field flow produced by a force F applied 
at the position (0, 0, z') can be expressed in the form 

u(r) = §7T- V 1 ^ -3 ^ - z ) z 'i H - z ')V||V||(-lnp) -F + o{e~ p/H ). (39) 
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The above relation can also be obtained by a direct expansion of the boundary- 
value problem in the small parameter H/ p j2?J . 

Relation (39) indicates that the correction to the far-field 0(p~ 2 ) asymptotic 
behavior of the fluid velocity u decays exponentially with p. Moreover, the 
vertical component of the force F does not contribute to the 0(p~ 2 ) behavior. 

Equation (39) can be rephrased as an expression for the asymptotic form T as 
of the full Green's function (6) 

T as (r,r') = -^-^H-hiH - z)z'{H - z')V\\V\{-ln\p - p'\), (40) 

where r' = p' + z'e' z . One of the gradient operators in the above formula has 
been applied to the primed coordinates to emphasize the Lorentz symmetry 
of the Green's tensor 

T as (r,r') = T ast (r',r) (41) 

(where the dagger denotes the transpose of the tensor). Due to the transla- 
tional invariance of the system in the lateral directions, the Green's function 
(40) satisfies the identity 

T as (r - q, r' - q) = T as (r , r) , (42) 

where the vector q has only lateral components. 

Using Eqs. (32) and (33) and noting that the Green's function (40) is quadratic 
both in primed and unprimed transverse variables, we find the relation 



T as (r,r') 



oo 



nr/H 3 



E' v£-(rX + *(rO, P>P\ 

(43a, b) 

oo 

E' V m + ( r KT*( r ')> P<P\ 



m=—oo 



which is analogous to the diagonal representation of the Oseen tensor (13). 
Equation (43a) combined with the displacement theorem for the Hele-Shaw 
basis fields (38) and identity (42) yields the symmetric representation of the 
asymptotic Green's tensor (40) 



oo oo 



T as (r, r') = £ ' £ ' v- + (r _ ei )S: y T(Q ij; m \ m')v™+*(v> - Qj) 



m=—oo m'=— oo 



(44) 
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4-4 Asymptotic form of the two-wall Green's matrix 



The asymptotic form of the matrix elements (26) can be obtained by projecting 
relation (44) onto the reciprocal basis fields w^ CT centered at the points Rj 
and Rj. The resulting expression involves the matrix elements 

(OwLrW I <' + «> = Smm'CiZ,; Ima), (45) 

where v^, + (i) = vj^, + (r — Qi). The elements (45) are diagonal in the azimuthal 
number m by cylindrical symmetry, they depend only on the vertical coordi- 
nate Zi of the point Rj = Qi + Zi& z , and they are real. Using these properties, 
the following asymptotic form of the wall Green's matrix (26) is obtained 

G^Alma | I'm'a') = —C{Z i -lma)S+\{Qi j -m\ m')C{Z^l'm! a'). 

nr]H 6 

(46) 

Due to the symmetric structure of the expression (46) and the symmetry 
property (37) of the scalar displacement matrix S^[, the Lorentz symmetry 

Gfjilma | I'm'a') = G™*{l'm'a' | Ima) (47) 

is manifest. We note that the presence of the Heaviside step function in relation 
(35) implies that 

Gf^lma | I'm'a') = for mm > 0. (48) 

The physical interpretation of the matrix C follows from the expression 

V- + (r - Qi) = E v w(r " Ri)C(Z i; Ima), (49) 

la 

which results from Eqs. (16) and (45). The matrix C(Z;lma) thus describes 
the transformation from the representation of the flow in terms of nonsingular 
Hele-Shaw basis v^ + (r — g { ) centered at the lateral position Qi to the spherical 
representation (12b) centered at R«. 



4 ■ 5 Multipolar flow fields 

An alternative interpretation of the matrix C is obtained by considering the 
far-field flow 

O - g 2 ; Z 2 ) = J T as (r, rXV - R 2 )w+ a (r' - R 2 ) dr' (50) 
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produced in the space between the walls by the multipolar force density 



F(r') = 5 a S (r'-R 2 )w+ mCT (r'-R 2 ) (51) 

centered at R 2 . By inserting representation (43 a) specified for the shifted 
asymptotic Green's function (42) with g = q 2 into (50) and using definition 
(45) of the matrix C we find that 

iC CT (r " Q2; Z 2 ) = ~^p<~(r ~ Q2)C(Z 2 ; Ima). (52) 

Thus, the matrix element C(Z 2 ; Ima) represents the amplitude of the Hele- 
Shaw basis field in the far-field multipolar velocity (50). Only one term 
contributes to this flow according to Eq. (52) because of the cylindrical sym- 
metry of the problem. 

The asymptotic multipolar flow fields (50) can also be expressed in terms of 
the matrix elements (46). To this end the right side of Eq. (50) is expanded in 
the spherical basis fields (12b) with the help of identity (16). The expansion 
yields the relation 

u^(r - q 2 ; Z 2 ) = J2 v+ mV ,(r - ROG-tf'mV | Ima), (53) 

I'm'a' 

where G^f is given by Eq. (26) with the Green's function T replaced with 
T as . The above expression relates the asymptotic flow uf^ la centered at the 
position R 2 and the spherical basis fields centered at a different position R x . 

We note that for each m only several force multipoles (51) produce a nonzero 
far- field velocity (50). This behavior results from the properties of the matrix 
C(Z 2 ; Ima) that appearers in relation (52); the form of this matrix is analyzed 
in Sec. 4.6 below. A further discussion of the multipolar fields in the space 
between the walls is given in Appendix B. 

4-6 Explicit expressions for the transformation matrix C 

A general structure of the matrix C can be inferred using scaling arguments. 
According to Eq. (12b) spherical basis fields v^ CT (r — Rj) are homogeneous 
functions of the order I + a — 1 of the relative-position vector r, = r — Rj. 
Similarly, Eqs. (31b) and (32) imply that the Hele-Shaw basis fields v^ + (r — 
Qi) are combinations of homogeneous functions of the order \m\ — 1, |m|, and 
\m\ + 1 of Tj. Since the coefficients C{Z i - 1 Ima) are independent of Tj, relation 
(49) implies that the non-zero elements of C[Zi \ Ima) satisfy the condition 

l + a-\m\<2. (54) 
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A detailed analysis of relation (49) reveals that the nonzero elements of the 
matrix C can be written in the form Uj 



C(Z;l±n o) = B±M Z), 



\m\ > 1, 



where B^ a (fx; Z) are the elements of the 3x3 matrix 



(55) 



-Z(H-Z) 



-Im(H-2Z) 
(/x + l)(2 /U + 3) 1 /2 



T{H - 2Z) 



±2/i 

(/i + l)(2^ + 3)V2 



2/x( M + l) 



1/2 



(/i + 2)(2/i + 3)(/i + 5) 1 /2 



with 



A ± ( A i) = ( T 2)" A i! 



4tt 



1/2 



(56) 



(57) 



(2// + l)(2/i)! 

The range A = 0, 1, 2 of the index A = I — \m\ in equation (56) result from the 
conditions \m\ < I and (54). All other elements of the matrix C vanish. 

We close our theoretical considerations with a remark that the asymptotic 
form Gfj of the Green's matrix Gij is approached exponentially for i?^- — > oo 
because Liron-Mochon formula (39) is exponentially accurate at large dis- 
tances. As in relation (39), the lengthscale for this approach is set by the wall 
separation H. The asymptotic expression (46) should thus be very accurate 
when the interparticle distance Rij is larger than several wall separations H. 
This conclusion is supported by our numerical results discussed in the follow- 
ing section. 



5 Numerical examples 



5.1 Matrix elements 



A typical behavior of the Green's matrix is illustrated in Figs. 1 and 2. 
The results are shown for the matrix elements 

Gi 2 (l -10 | l/itr) = (-l) a G* 2 (110 | I -na) (58) 

versus the lateral distance Q12 scaled by the wall separation H for several values 
of the parameters /, cr, and ji > 0. The elements (58) play a special role in our 
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theory because in the asymptotic regime Qi 2 ^ H they are directly related to 
the far-field multipolar flow (50) according to Eq. (B.7). In all examples, the 
vertical coordinates of the points (1) and (2) are Z\ = Z 2 = \H. For other 
configurations the matrix elements have a similar behavior. 

We present our results in the form of the rescaled elements defined by the 
relation 

G 12 (l =f10 | I ±fia)=H»- l -° +1 G 12 (l =pl | I ±p)$± ((J+1) (ei 2 ). (59) 

For those values of parameters / and \i for which the matrix elements (59) 
do not vanish, the factor $±( M+1 )(£i 2 ) ~ Q12 corresponds to the far-field 
behavior of Grf|(l =FlO | / ±/xcr), according to Eqs. (35) and (46). In the 
asymptotic regime Qu H the rescaled elements G as i2(l =Fl | / ±fia) de- 
pend only on the vertical coordinates Z\ and Z 2 . The nonzero asymptotic 
elements are quadratic functions of the vertical coordinate Z\, and they are 
at most quadratic in Z 2 but can also be linear or constant in this variable, 
as indicated by Eqs. (46) and (56). The far-field flow (50) is related to these 
elements by Eq. (B.9). 

Figure 1 illustrates the behavior of the matrix elements 612(1 —10 | I la) with 

I + a < 3. (60) 

All these functions approach nonzero asymptotic values G as i2(l —10 | 11 a) 7^ 
for large interparticle distances Qu 3> H, according to Eqs. (46) and (54). 
The corresponding behavior of the unsealed matrix elements (58) 

G 12 (l =pl I I ±1 a) ~ q u 2 , g 12 > H, (61) 

follows from Eqs. (31a), (35), and (46). 

The matrix elements (61) are directly related to the multipolar flow fields (50), 
as indicated by Eq. (B.7). Therefore, we find that Eq. (61) corresponds to the 
slowest possible far-field decay of the flow produced by a multipolar force 
distribution. The multipoles (60) include the horizontal Stokeslet (I — 1, a — 
0), rotlet (I — 1, a — 1), stresslet (I — 2, a — 0), and three other multipoles, 
one of which has the spherical-harmonics order 1 = 3. The numerical results 
shown in Fig. 1 indicate that the approach of G± 2 to the asymptotic values is 
exponential, which is consistent with our discussion in Sec. 4. 

Figure 2 illustrates the behavior of matrix elements (58) for I — 4 and a = 0. 
Unlike the elements presented in Fig. 1, the rescaled elements Gi 2 (l — 1 | 4 m 0) 
with m — 0, 1 vanish at large separations g ^> H, which is consistent with 
condition (54). The remaining rescaled elements with m = 2,3,4 tend expo- 
nentially to nonzero asymptotic values. 
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1 = 2, (7 = 1 



| = 3, cr = 



12 3 

Pl 2 /i? 



Fig. 1. Rescaled matrix elements Gi2(l —10 | I la) versus lateral distance Q12 for 
Z\ = Z2 = \H; values of parameters I and a as labeled. Exact solution (solid lines); 
asymptotic behavior (dotted lines). 

5.2 Applications in multiparticle hydrodynamic-interactions algorithms 



The simplest numerical application of our asymptotic formulas (46), (55), and 
(56) is to implement them directly in the induced-force-multipole equation 
(21). To this end, the matrix (26) is represented as the superposition of the 
long-range asymptotic part and the short-range correction, i.e., 

Gijilma I I'm' a') = Gf^lma \ I'm' a') + bG^ilma \ I'm' a'). (62) 

The asymptotic part G^(lma | I'm' a') can be evaluated from our explicit for- 
mulas at a low numerical cost. To obtain the correction term 5Gij(lmo~ \ Ima), 
first the expression Gij(lmo~ \ I'm' a') is calculated using the Cartesian- repre- 
sentation method described in 

Urn 

and next, the asymptotic expression 
is subtracted from the result. Since the correction term is short ranged, the 
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Fig. 2. Rescaled matrix elements Gi2(l —10 | 4m0) versus lateral distance Q12 
for Z\ = Z2 = \H; values of parameter m as labeled. Exact solution (solid lines); 
asymptotic behavior (dotted lines). 



matrix 5Gij(lmcr | I'm 1 a 1 ) can be truncated by setting 

5Gij{lma \ I'm! a') = for > £> as , 



(63) 



where the truncation distance g as is of the order of several wall separations H. 
The results shown in Figs. 1 and 2 and other similar tests indicate that the 
asymptotic approximation for the Green's matrix is very accurate for g as > 
3H. Thus, the numerically expensive contribution 5Gij has to be evaluated 
only for the neighboring particles at an O(N) cost. 

To test our asymptotic approach and illustrate the role of the long- and short- 
range contributions to the Green's matrix G, we consider a benchmark case 
of a linear rigid array of N touching spheres translating in the center plane 
in the space between closely spaced walls. The spheres are on a line parallel 
to the x direction and the array is moving either in the x (longitudinal) or 
y (transverse) direction. We focus on the translational friction coefficients 
evaluated per particle, 



act 
C 



N 

1 E 



/•tt act 



a 



(64) 



tt 



where Qj aa is the aa component of the translational resistance tensor C, 
defined in Eq. (2), and G is the one-particle lateral translational resistance 



coefficient lla 14 



As illustrated in Fig. 3 (see also discussion in [3, 3]) the longitudinal and 
transverse friction coefficients (64) behave differently. The longitudinal co- 
efficient (q x decreases with the length of the array N, while the transverse 
coefficient (j^ increases with N. For tight configurations with small gaps be- 
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Fig. 3. Normalized longitudinal and transverse friction coefficients per particle £q x 
and CJc for linear arrays of touching spheres in the center plane between the walls, 
for wall separation H/2a = 1.05. Crosses represent exact results; open symbols 
represent asymptotic approximation (63) with Qas/H = 1 (squares) and 2 (circles); 
solid symbols correspond to truncation (66) with Qq/H = 2 (squares), 4 (circles) 
and 6 (triangles). For the longitudinal coefficient (q x the open squares and circles 
coincide with the crosses. 

tween the wall and particle surfaces (the case show in Fig. 3) the decrease of 
(q x is moderate because the friction force is dominated by the local resistance 
due to the dissipation in these gaps. In contrast, the increase of (q 1 is large 
due to collective long-range effects. 

The mechanism of these collective effects can be explained using the results 
for the pressure field around arrays of the length N = 10 and 20 plotted in 
Figs. 4 and 5. The figures show the normalized asymptotic far-field pressure 

pas = HiriU)- 1 ^ (65) 

(where U is the velocity of the array), which was evaluated using the method 
described in Appendix C. The results for the longitudinal motion of the array 
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Fig. 4. Contour plots of the normalized asymptotic pressure field (65) around a 
rigid array of N = 10 touching spheres moving in the center plane between parallel 
walls, for wall separation H/2a = 1.05. Longitudinal motion (a); transverse motion 
(b). The dotted lines delimit the regions where the asymptotic approximation is not 
valid. 

shown in Figs. 4(a) and 5(a) indicate that the pressure field is only weakly 
affected by the length of the array, and its magnitude is the largest near 
the array ends. In contrast, the pressure shown in Figs. 4(6) and 5(6) for 
the transverse motion increases approximately linearly with the array length 
N, and its magnitude is maximal near the chain center. This large pressure 
amplitude is associated with the flow of the displaced fluid around the ends of 
the array in the confined space. The flow is significant over the distance that 
scales with the length of the array / = 2Na (where a is the sphere radius). 
In the Hele-Shaw regime the pressure gradient is proportional to the fluid 
velocity; hence, the pressure itself is proportional to N. 

To further elucidate the effects of the short-range and far-field flow compo- 
nents, the exact numerical results for the resistance coefficients (64) are com- 
pared in Fig. 3 with the asymptotic approximation (62) and (63). We also 
show results obtained using a much cruder approximation, where the whole 
Green's matrix is truncated at a certain distance p , i.e., 

Gijilma | I'm a') = for ^ > g . (66) 

Our numerical calculations indicate that the truncation (66) yields poor re- 
sults. The far-field flow contribution is especially important for the transverse 
motion of the array because of the positive-feedback effect: For this motion 
the dipolar Hele-Shaw flow field generated by a given particle acts as a back 
flow on the other particles. This back flow, in turn, produces an increase in 
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Fig. 5. Same as Fig. 4, except that for N = 20. 

the induced force distribution that generates the dipolar flow. This back flow 
mechanism, resulting in the large transverse resistance, is consistent with our 
discussion of the pressure field shown in Figs. 4(b) and 5(b). 

In contrast to the crude approximation (66), a truncation of the short-range 
part (63) of the Green's matrix yields accurate results already with moderate 
values of the truncation parameter g as . The results shown in Fig. 3 indicate 
that the truncations at g as /H = 1 for the longitudinal motion and at £> as / H = 
2 for the transverse motion are sufficient. The results with g as /H > 3 (not 
shown) are essentially indistinguishable from the exact results. 



6 Conclusions 



Our paper presents a complete analysis of the far-field flow produced by an 
arbitrary force multipole in the space bounded by two parallel planar walls. 
We have shown that a force multipole characterized by the multipolar num- 
bers Ima produces, at large distances, a Hele-Shaw flow driven by a two- 
dimensional multipolar pressure field of the azimuthal order m. The amplitude 
of this flow has been explicitly obtained for an arbitrary order of the source 
force multipole. 

Our asymptotic results were applied to evaluate the multipolar matrix el- 
ements Gij(lma | I'm' a') of the Green's tensor for Stokes flow in the wall- 
bounded domain. This matrix is used in our recently developed algorithm 
[Til 15| for evaluation of the multiparticle friction tensor £y in a suspension 
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confined between two parallel walls. The elements of the matrix Gij are equiv- 
alent to the expansion coefficients in the displacement theorem for Stokes flow 
in the bounded domain. Such a displacement theorem relates the flow pro- 
duced by a force multipole centered at a point Rj to nonsingular multipolar 
flows centered at a point Rj. We have shown that in the far- field regime the 
matrix elements Gij(lma \ I'm' a') can be expressed in terms of much simpler 
displacement formulas for the two-dimensional scalar potential. 



We have found that the matrix achieves its asymptotic behavior when 
the lateral distance between the centers of the particles i and j exceeds sev- 
eral wall separations H. Evaluation of the exact matrix elements in terms of 
lateral Fourier integrals derived in 0, 14 1 is thus needed only for the neigh- 
boring particles. Therefore, application of the asymptotic expressions in our 
hydrodynamic-interaction algorithm yields an important improvement of its 
numerical efficiency. (The far-field contribution to the Green's matrix cannot 
be simply neglected — for some problems such a crude approximation leads to 
entirely wrong values of the friction matrix; cf. discussion in section 5.2). 



Several other important consequences stem from the fact that we have reduced 
a complex hydrodynamic problem to a simpler problem of a two-dimensional 
scalar potential. First, since for a scalar potential the multipolar flow fields in 
a periodic system are known |30[|, the results of our analysis can be used to de- 
velop an algorithm for hydrodynamic interactions in a periodic wall-bounded 
system. Without the asymptotic expressions, evaluation of the periodic hydro- 
dynamic Green's matrix would be much more difficult, as discussed in [27]. 



Next, for scalar potentials, fast multipole and PPPM acceleration techniques 
are well developed [3l(. Combined with our asymptotic results, such meth- 
ods can be applied for fast evaluation of hydrodynamic interactions in wall- 
bounded suspensions. Development of accelerated algorithms for suspensions 
will require implementation of certain techniques that are specific to mul- 
tiparticle hydrodynamic systems, e.g. an appropriate preconditioning of the 
Green's matrix and incorporating the lubrication interactions into the calcula- 
tion scheme. These techniques were used in accelerated Stokesian-dynamics al- 
gorithms for unbounded suspensions [32| and [33| . Our present asymptotic re- 
sults greatly facilitate development of accelerated algorithms for wall-bounded 
systems, and our research is currently focused on this problem. 



S. B. would like to acknowledge the support by NSF grant CTS-0201131. E. W. 
was supported by NASA grant NAG3-2704 and, in part, by KBN grant No. 
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A Spherical basis 



The spherical basis sets of Stokes flows vf m(7 used in the present paper are 
normalized differently than the corresponding sets v^^ FS ^ introduced in |2ll |. 
The transformations between the basis fields are 



v w( r ) = N ia n im v , 



1^-1, -(CFS), 



Ima 



(A. la) 



where 



and 



vL,(r) = AW>:T S) (r) 



Nio = 1, 



N n = -(l + l)-\ 
N l2 = l[(l + l)(2l + l)(2l + 3)]-\ 



ni % 



47r (l + m)\ 
21 + 1 (l-m)\ 



1/2 



(A.lb) 
(A.2a) 
(A.2b) 
(A.2c) 

(A.3) 



Below we list the explicit expressions for the angular coefficients 
for spherical basis fields (12) in our present normalization, 



V 



ImO 



(21 + If 



l(2l-lf l ll ~ lm ~2 



(A.4a) 



and 



V 



V 



Iml 



-liY, 



l + l 

^lm2 — AYn+lm 



llmi 



V 



ImO 



atiY, 



i l—lmi 



V 



Iml 



l + l 



JiY, 



llmi 



I 



lm2 



2(2/ + 1) 



otiYn-im + 



I 



(Z + l)(2/ + l)(2Z + 3) 



AY 



I l+l mi 



where 



Y„_ lm (f) = a7 l r- l+l V 



Y 8+lm (r) = /3fV +2 V 



r < l+1 ^Y lm ir: 



Y Hm (r) = 7; r X V s Y" Zm (r) 



(A.4b) 
(A.4c) 
(A.5a) 

(A.5b) 

(A.5c) 

(A.6a) 
(A.6b) 
(A.6c) 
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are the normalized vector spherical harmonics, as defined by |34|. Here 

IW*) = n^(-l) m Pr{coB0)^ (A.7) 
are the normalized scalar spherical harmonics, and 

ai = [1(21 + l)] 1 / 2 , (A.8a) 

A= [(/+1)(2/ + 1)] 1/2 , (A.8b) 
7i = -i[i(i + l)] 1/8 . (A.8c) 



B Flow fields generated by force multipoles 

In this Appendix we express the flow field 

u lma (* - g 2 ; Z 2 ) = J T(r, r')5 a S (r' - R 2 )w+ CT (r' - R 2 ) dr' (B.l) 

produced in the space between the walls by the force multipole (51) in terms 
of the elements of the Green's matrix G\ 2 . Using relation (16) to expand the 
right side of Eq. (B.l) into the non-singular spherical basis fields (12b) yields 

u /mCT (r - q 2 - Z 2 ) = ]T v+ mV ,(r - R 1 )G 12 (/'mV | Ima) (B.2) 

I'm'a' 

where the definition (26) of the Green's matrix have been used. Equation (53) 
is the asymptotic form of the above relation. 

Equation (B.2) can be simplified by evaluating it for r = R x and noting that 

v+ m , ff ,(0) = for l' + a f >l (B.3) 

and 

vL'o(0) = (|vr)- 1/2 e m ,, (B.4) 

where m! = 0, ±1, and 

e±i = =F-^(e x . ± ie y ), e = e z , (B.5) 

which follows from Eqs. (A. 5). According to the above expressions, only three 
terms of the sum (B.2) contribute to the result, 

iw(Ri - Q 2 ; Z 2 ) = (fvr)- 1 / 2 Y, G 12 (lm'0\lma)e m ,. (B.6) 

m'=-l,0,l 
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In the asymptotic regime Q\ 2 ^> H this relation simplifies because the asymp- 
totic matrix elements (46) vanish for m'm > 0, as discussed in Sec. 4. Taking 
this observation into account we thus find 

u?V(Ri - g 2 ; Z 2 ) = (f vrr 1 / 2 ^! =pl | / ±/i a)e Tl (B.7) 

where \x = 1,2,... The dependence of the flow field (B.7) on the lateral 
relative-position vector g 12 can be made more explicit by using relation (59) 
and the identity 

$ ±(M+i)(£i2)e T i = T2' 1/2 (fi + lr'V^teia), (B.8) 

which yields 

ur ±M(7 (Ri-f? 2 ;Z 2 ) = T(|7r)- 1 / 2 (//+l)- 1 //^- CT+1 G- 12 (l =pl | I ±»a)V^(g 12 ). 

(B.9) 

We note that relation (B.9) is consistent with (52) due to Eq. (32). 



C Far-field pressure distribution 

As discussed in Sec. 4, the flow and the pressure fields in the Hele-Shaw 
asymptotic regime (27) are uniquely related (up to an additive constant in the 
pressure). Thus, many of the asymptotic formulas, expressed here in terms of 
the velocity fields, can be translated into the corresponding expressions for 
the pressure. 

This remark applies, in particular, to Eq. (52) for the asymptotic multipolar 
flow (50). We introduce the asymptotic multipolar pressure field p^(r), which 
is defined by the relation 

<Ur - f? 2 ; Z 2 ) = -\rf x z{E - z)V\\pT m(7 {P ~ fel Z 2 ). (C.l) 

Using Eqs. (30) and (32), the flow-field identity (52) can be transformed into 
the corresponding pressure identity of the form 

pZ a (P - Q2; Z 2 ) = -Jp® m (P ~ Q2)C(Z 2 ; Ima). (C.2) 

Equation (C.2) can be conveniently used to evaluate the far-field disturbance 
pressure p as in a many-particle system. This equation describes the asymp- 
totic pressure produced in the far-field regime by a single force multipole 
(51), as indicated by Eqs. (50) and (C.l). To determine p as , the multipolar 
moments fi(lma) of the force distributions (18) induced on the surfaces of 
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particles i = 1, . . . , N are evaluated by solving the force- multipole equations 
(21). Combining the solution with (C.2) yields 

N 

f(p) = EE'ftW^-ft). (c.3) 

i=l m 

where 

Qi(m) = -^E^^iH/.W. (C4) 

la 

The contour plots in Figs. 4 and 5 were obtained using this method. 
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